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Abstract: Capturing disused satellites in orbit and their controlled reentry is the aim of the 

DEOS space mission. Satellites that ran out of fuel or got damaged pose a threat to working 

projects in orbit. Additionally, the reentry of such objects endangers the population as the place 

^TJ^ of impact cannot be controlled anymore. This paper demonstrates the modelling of a rendezvous 

szenario between a controlled service satellite and an uncontrolled target. The situation is 

I— I modelled via first order ordinary differental equations where a stable target is considered. In 

Q_J order to prevent a collision of the two spacecrafts and to ensure both satellites are docked at 

/^ the end of the maneuver, additional state constraints, box contraints for the control and a time 

dependent rendezvous condition for the final time are added. The problem is formulated as an 

r^ optimal control problem with Bolza type cost functional and solved using a full discretization 

5^ approach in AMPL/IpOpt. Last, simulation results for capturing a tumbling satellite are given. 
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Kj^ 1. INTRODUCTION docking and berthing with a noncooperating target and 

fvq the controlled reentry into atmosphere in coupled condi- 

QQ tion, see Scllmaier et al. (2010). 

1^ According to NASA's Orbital Debris Program, see NASA This paper presents a model to calculate the optimal con- 

\Q Orbital Space Debris Program (2011), approximately trol schemes for a DEOS capturing mission. In particular, 

• 19.000 known objects with diameters larger than ten cen- the equations for a rendezvous maneuver to an uncontrol- 

y timeters and estimated 500.000 objects between one and lable target are derived and solved for one example with a 

^ ten centimeters are currently in orbit. Due to their average direct approach. The presented results show the feasibility 

,__! speed of 10 km/s even such small objects may cause serious of such a mission. Additionally, reference trajectories and 

•• damage to working satellites or the international space sta- control schemes for different situations are obtained which 

^ >^ tion. Additionally, with each collision the number of debris are extendable to the case of a tumbling target. 

k> increases rapidly as the cloud of around 600 fractions of The following section deals with the derivation of the 

^^ the 'Cosmos 2251' and 'Iridium' crash in 2009, see Pickup optimal control problem. The equations of motion for the 

C^ (2009). relative position and orientation, the terminal condition. 

Aside from orbital collisions, an uncontrolled reentry into the state and control contraints and the cost functional 

the atmosphere poses a severe threat. As the entrance are presented. Afterwards a direct approach is applied to 

corridor is highly difficult to precalculate, with every crash goive the optimal control problem. In the fourth part the 

pieces of beryllium or titan may reach the earth's surface numerical results for a flyaround maneuver with a stable 

and hit populated areas. In October 2011 the german spinning target are presented and discussed. Thereafter 

satellite 'Rosat' crashed uncontrolled and even a few hours a short introduction in simulating the trajectories in a 

prior to impact experts were unsure where fractions will virtual reality is given. An outlook of possible extensions 

strike the earth, see KHnkrad (2011). ^ ^ ^ concludes the paper. 

The planned german technology demonstration mission 

DEOS (Deutsche Orbitale Servicing Mission) takes the 2. SETUP OF THE RENDEZVOUS MODEL 
next step towards capturing and controlling disused satel- 
lites. Considered as a feasibility study, DEOS consists of The motion of each spacecraft consists of two independent 
two spacecrafts, one acting as uncontrollable target and 3D subsystems, the relative position and the relative ori- 
one as servicer. The aims of the mission are to demonstrate entation. The equations of motion for the relative position 



are the so called Clohessy-Wilshire equations. Further 
information on orbit analysis and design can be found 
in Schilling (2011). Note that throughout this work we 
shorten notation by omitting the time dependency if no 
specific time instant is considered. 

2.1 Relative position 

Assuming circular orbits the subsystem of the relative 
position is derived via Kepler's equation for the two 
body problem and a linear taylor expansion, see, e.g., 
Alfriend et al. (2009). As a result, a second order ordinary 
differential equation system for the position of a body 
relative to a moving reference point in orbit is obtained. 
Throughout this work, we consider the uncontrolled target 
to be that reference point. Apart from the standard state 
variables x, y and z denoting the relative distance in 
radial, direction of movement and out of orbit component 
respectively, we introduce the additional state variables Vx, 
Vy and Vz denoting the respective velocities to transform 
the system into six first order differential equations. 

X = Vx (1) 

y = vy (2) 

z = Vz (3) 

Ux 



2nVy + 3n a; + 



Vy = —2nVx 



-n z + 
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m 



(4) 
(5) 
(6) 



Here m denotes the mass of the servicer spaecraft and n 
is the mean motion of the reference point, i.e. the angular 
velocity in orbit. For the dynamics it can be seen that the 
z coordinate is decoupled from the x-y-system, i.e. the out 
of orbit motion is independent of the motion in the orbit 
plane. If no controls Ux, Uy, Uz are applied an analytical 
solution is known, again see Alfriend et al. (2009). Yet for 
non vanishing controls no analytical solution is known. For 
this reason, we evaluate the dynamics numerically. 

2.2 Orientation 

In space mission modelling one commonly uses the quater- 
nion representation to derive the equations for the orienta- 
tion of a satellite. Quaternions are a special parametriza- 
tion of the Euler-axis/angle description which is nonsin- 
gular and continuous. A quaternion can be considered as 
a four dimensional vector q :— [qi, q2, qs, q4]^ . The first 
three components are called the vector part and the fourth 
is the scalar component. The denomination results from 
the fact that a quaternion represents a rotation around an 
axis in three dimensional space, i.e. the vector part, by 
a specific angle given by the scalar component. To get a 
unique description of the orientation, the quaternions are 
supposed to be of length one, so called unit quaternions. 
Further information on quaternions and their relation to 
other rotation types can be found in Tewari (2007). 
The first derivative of quaternions is given by 

1 ' 



Q = 



q, 



where lo 



[lUx, LUy, uJz]^ is the vector of the angular 
velocities and (g) represents the quaternion multiplication. 
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see also Stevens and Lewis (2003) for a derivation. In its 
components the dynamics of the quaternions are of the 
form 



ae{^,r}, (7) 



where S and T denote the service and target spacecraft 

respectively. 

Note that the angular velocity of a tumbling object is 

typically not constant with respect to time. Denoting the 

inertia tensor with J we incorporate that aspect using 

Euler's gyroscopic equation 

J-w-|-J-a;-|-a;x (J-u;)= m, 

which provides the needed information, see, e.g., Chobotov 
(1991). Under the assumptions that mass distribution is 
constant over time, i.e. J = 0, and the regarded system 
coincides with the principal axis of the body, that is 
Jik = for all 
equations 



=/= k, we obtain the set of first order 
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(8) 
(9) 

(10) 

(11) 
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where m := [mx^ nriy, ruz] defines the momentum 

vector. Note that no momentum m is applied on the right 

hand side of the target satellite as it is supposed to be 

uncontrolled. 

With these equations of motion we obtain a system of 

twenty first order ordinary differential equations (1)~(13). 

Here, we denote the combined state vector by 

X := [x, y, z, Vx, Vy, Vz, wf , u^ , wf , wj, wj, 

'^f. ^f, qi, Qii (1a, ql, ql, qI, qIV 

and the control vector by 

- r iT 

u . \vix., n^i, u^. rrix, TTir,,. Tfi^i 



2.3 Terminal condition 

To model the docking of the two spacecraft some terminal 
constraints have to be introduced. We define a docking 
point for each satellite in their body fixed coordinate 
system, dP := [df, df , df]^ and (f := [d^, d^, df']^, 
whose local positions and velocities have to match at the 
end of the maneuver. Similar conditions can be found in, 
e.g., Alfriend et al. (2009) and Boyarko (2010). To get a 
common representation, all components are converted to 
the earth centered inertial coordinate system. We denote 
the corresponding transformation matrix from the rotated 
body fixed coordinates to the inertial system by R and 
add the index £ to variables defined in the earth central 
inertial coordinate system. Computing the transformation 



matrix from the components of the related quaternion, see, 
e.g., Tewari (2007), we obtain 
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(14) 

Multiphcation with this matrix performs the transforma- 
tion of a vector from the unrotated coordinate system 
to the rotated one. The inverse transformation, as we 
need it in this case, can be done by muhiphcation with 
the transposed matrix of R. Applying the transformation 
for both the service and the target satellites, the local 
difference of the docking point is given by 



4-4 = r 



R 



'^d^ - R^ 



^d\ 



(15) 

with r :~ [x, y, z]^ being the cartesian distance of the two 
spacecrafts. Similarly, the angular velocity of each satellite 
in the earth centered inertial frame is given by 





= W 



a e {S, T} . 



Hence, we obtain the difference of the velocity via 



(i?^^d^) - u;J X (if^d^ 



(16) 

To simplify this equation we assume that the angular 
velocity and the orientation of both spacecrafts coincide 
at the end of the time interval, i.e., 

q^(tf)~q^{tf)^Q (17) 



\h)-^''^^i) 



0. 



(18) 



Note that the docking point definition is required to be 
conform with these assumptions. To illustrate this point, 
consider the docking points to be located at the back 
of the satellites. Then, since the satellites would have to 
overlap, there exists no physically possible configuration 
where these points and the satellite orientations coincide. 
Assuming (17) and (18) to hold we obtain 

and the docking conditions (15) and (16) simplify to 
R{tf)'^ {(f -d^\-r = 
u>s{tf) X R{tf)^{d^ -d^)-r = Q. 

2.4 State and control constraints 

As pointed out in the previous subsection, the case of 
colliding satellites has to be treated. To exclude such an 
occurance, we first introduce a spherical safety area around 
each satellite with radius r^' > and r^ > respectively. 
In our case, the usage of spheres is reasonable as the 
geometry of the satellites is supposed to be simple. Based 
on the safety areas, we add the state constraint 

x^ + y^ + z^ > (r^ + r'^f 

to the optimal control problem which resolves the possible 
collision issue. Note that the safety area needs to be 
chosen carefully since the state constraint together with 
the terminal condition may exclude the existence of a 
solution. 



Apart from the states we also consider the controls to 
be constrained. In particular, we suppose that the thrust 
vector u :— [u^, Uy, Uz]^ is bounded via 

ul+Uy+ul< Umax 

and each momentum is bounded via 



mj < rrin 



X, y, z. 



The choice of the limitations in the momentum control is 
motivated by the fact that the satellites attitude control 
system can be realized via control moment gyroscopes. 
These can cause a maximal momentum nimax around each 
axis. The bound of the thrust vector is due to future 
extensions of the model where only one thruster is installed 
and the spacecraft has to rotate to accelerate in a specific 
direction. 

As the thrusters are mounted on the satellite their direc- 
tion depends on the current attitude of the spacecraft. 
Therefore we transform the calculated thrust u^^ Uy and 
Uz from earth centered coordinates into the body fixed 
system. The transformation is done by multiplication with 
the rotation matrix R defined in (14) and reveals 
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where mi, U2 and u^ denote the controls in the body fixed 
coordinates. Note that as R is orthogonal the length of u 
remains unchanged by this transformation, and hence the 
upper bound ■Umax on the controls in the body fixed and 
the earth centered coordinates are identical. 

2.5 Cost functional 

Throughout this work we consider minimizing the follow- 
ing Bolza type cost functional over the set of allowed 
control functions ■u(i) :— [u{t), m{t)]^ and terminal times 
tf. We define the cost functional as follows 

J{u{t),tf) = Itjtf + r lu \\u{t)\\l + l„, \\m{t)\\ldt. 

(19) 
This functional represents a combination of the terminal 
time, the position controls and the orientation controls. 
The non-negative constants It, , lu and l-m can be used to 
force scenarios such as time-optimalty or minimal control 
effort. Note that the 2-norm || • II2 has to be used if results 
shall hold for both controls (ui, U2, u^)^ and (uj., Uy, w^)^. 

3. DISCRETIZATION AND IMPLEMENTATION 

To solve the optimal control problem 
Minimize J{u{t),tf) 
subject to the dynamics (1) — (13) 
with initial and terminal conditions 

a;(0) = xo 

= q^{tf)-qS{ts) 
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ulit) Vte[0,t/] 



(20) 



m,„ax > |m-i(i)| i^x,y,z ViG[0,t/] 

x^{t) + y^{t) + z\t)>{r'^ + r^f Vie[0,t/] (21) 

a direct approach can be applied. Using a full discretiza- 
tion with fixed step size At — tf/N where N is the number 
of steps, the discretized dynamics of the state x and the 
discretized constraints (20), (21) form the constraints of 
the resulting finite nonlinear program. The optimization 
variable now contains not only the discretization of the 
control u and the free terminal time i/, but also the 
discretization of the entire state trajectory. With the ex- 
ponent k we denote the state of a variable at the fcth 
discretization point, i.e. x'' — x{k ■ At). 
Since the gyroscopic equation is stiff the implicit trape- 
zoidal method 
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fc = 0,l,. 
is used to discretize the dynamics of the state. 



,iV-l 

Note that 

the components of the angular velocity would soar even 
for very small step-sizes if, e.g., Euler's method was used. 
We implemented the discretized optimal control problem 
in the modelling language AMPL and solved it using the 
interior point optimizer IpOpt. For details on AMPL and 
IpOpt see, e.g., Fourer et al. (2003) and Wachter and 
Bicglcr (2006). For the discretized problem, the constraints 
are defined pointwise at the discretization points. To 
avoid constraint violations at intermediate time instances, 
adaptation techniques for the discretization grid may be 
used, see, e.g., Betts and Huffman (1998), or the bounds 
may be tightened. 
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Geometrically, the above conditions force the servicer to 
fulfill three tasks: It has to fly around the target, turn 
around to align the docking points and adopt the rotation 
such that the angular velocities coincide. 

Within the cost functional (19) we consider the weights 
/f, = 1, Z„ = 1 and Im = 1, i-C, the terminal time 
and the control costs are weighted equally. The remaining 
constants of the optimal control problem are displayed in 
Table 1. 



variable 



value 



description 



orbit radius [m] 

gravitational constant [N{m/kg)'^] 
mean motion [1/s] 
satellite mass [kg] 
"max U.it) maximum thrust [N] 

maximum torgue [Nm] 
Targ. angular mass around x [kg/m'^] 
Targ. angular mass around y [kg/m'^] 
Targ. angular mass around z [fcg/m^] 
Serv. angular mass around x [kg/rri^] 
Serv. angular mass around y [fcg/rra^] 
Serv. angular mass around z [kg/m'^] 
safety-area around target [m] 
safety-area around servicer [m] 
docking point target [m] 
docking point servicer [m] 

Table 1. Constants of optimal control problem 
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4. NUMERICAL RESULTS 

To illustrate our results we consider a flyaround maneuver, 
that is a situation where the service satellite is located 
on the wrong side of the target pointing with its docking 
interface towards it. Here, we assume that the docking 
interface of the servicer is located on the 'front' whereas 
the docking point of the target is at its 'back', e.g. the 
targets thruster where some kind of hook clasps. The 
radius of the safety area for each satellite is chosen to 
be of size one. For both satellites we assume that they 
are axially symmetric with respect to the y axis which is 
justified for the regarded satellite geometries. Note that 
the choice of the docking points enables us to use the 
simplified equations for the ternfinal condition (17) and 
(18) and creates no impossible final states. 

Upon start of the maneuver, we suppose that the servicer 
is in a non-rotated initial state and its angular velocity is 
zero. The target also starts in a non-rotating initial state 
but in contrast to the servicer it is rotating with a constant 
angular velocity of 3°/s around its y axis resulting in a 
stable motion with respect to time. The initial values are 
as follows: 
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To convert the optimal control problem into a finite 
nonlinear program we used a discretization with N = 370 
steps. The results of our calculations are displayed in 
Figures 1 and 2 which show the states and controls over 
time. For the optimal terminal time we obtain tf = 369.61 
seconds, which is about 6.16 minutes. 
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Fig. 1. Position and thruster control of the servicer 

From Figure 1 one observes that the position control 
never reaches the maximum. The fact that the allowed 
thrust is never applied is due to the chosen weights in 
the cost functional. The controls are smooth, except one 
point at about 190 seconds, but there occurs no bang-bang 
behaviour. Because of the state constraints the satellite 
cannot fiy to its terminal position in a straight line. In the 
beginning it accelerates in positive x direction - which is 
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constraint the flight path is much more direct and the 
servicer would crash into the target. 




Fig. 2. Orientation, angular velocity and momentum con- 
trols of target and servicer (dashed) 

due to the initial orientation - performed by a negative 
control of ui. Note that the applied control has to be 
interpreted with respect to the current attitude of the 
satellite, cf. Section 2.4. Upon termination of the maneuver 
the negative control in y direction represents the breaking 
to prevent collision with the target and to obtain zero 
relative velocity. As expected, the position in x and y 
nearly produce a symmetrical arc with respect to half the 
maneuver time. 

The orientations in Figure 2 are plotted in Euler-j/xz- 
angles with the angle notation (/), 9 and ^. One observes 
that the spin of the target around the y-axis is periodical. 
In the first part of the maneuver the main momentum 
control of the servicer is applied to perform the turning. To 
this end, the almost maximal momentum is applied around 
the z-axis and some additional correction in x direction. 
This control first accelerates the angular velocity of the 
servicer around the turning axis and then slows it down 
such that upon termination of the maneuver the angular 
velocity around the x- and z-axis are zero. After about 50 
seconds the momentum around the j/-axis is increasing. It 
reaches the maximal possible amount and nearly keeps this 
value until the end of the maneuver. This is the control for 
adapting the angular velocity around y. The acceleration 
of LOy is timed such that upon termination both the angular 
velocity and the orientation of servicer and target are 
equal. 

Computing the total control costs of the thrust and the 
momentum control, we obtain 



N-l 



"total 



Wtotal 



= At ^ {ulf + {ulf + [ulf = 15.7662 



fc=0 
JV-l 



= At ^ [mlf + [mlf + {m'lf ^ 295.5767. 
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Hence, the costs of the momentum control outweigh the 
costs of the position control for this particular example. 
The optimal value of the cost functional sums up to 

J = i/ + utotai + rritotax = 680.9548. 

To demonstrate that the calculated control depends on 
the state constraints we calculated the optimal trajectory 
with the same initial values omitting this constraint. In 
Figure 3 one can see that with the safety areas the servicer 
flies a nearly circular arc around the target. Without the 




Fig. 3. Trajectory of the servicer with state constraint 
(blue) and without state constraint (red) 

5. SIMULATION 

Since the computed results are hard to imagine from the 
plots, we developed a virtual reality where the calculated 
data can be imported and rendered images of the ma- 
neuver can be created. These images are then combined 
to a video showing the entire rendezvous maneuver. The 
virtual reality was created with the open source software 
Blender 3D. Using an additional import tool the trajectory 
data is bound to the corresponding satellite. Moreover, 
the calculated quaternions are used directly to rotate the 
models. For information on Blender see, e.g., van Gumster 
(2009). The satellite models were designed based on the 
first design studies of the german aerospace center, see 
Sellmaier et al. (2010). In Figure 4 one can see an image 
of the two satellites in the created virtual reality. 




Fig. 4. Rendered picture of the satellites 
6. OUTLOOK 

The shown example only covers a stable rotating target 
satellite, i.e. all angular velocity components are constant 
with respect to time. Yet, the presented model can also be 
used to calculate optimal solutions for a tumbling target 
which may occur if a satellite runs out of fuel or energy 
to keep its orientation. This case is under investigation 
at the moment to handle the additional difficulty of a 
moving target point. Note that it is possible to calculate 
the maximal tumbling motion which still allows a docking 
maneuver to be performed. If this maximum is exceeded, 
the service satellite can be fitted with an additional manip- 
ulator arm to grasp the target to establish a connection. 



To optimize such a so called berthing maneuver the model Wachter, A. and Biegler, L. (2006). On the implementa- 



derived in this work has to be modified by adding a model 
of a robot arm, see also Fehse (2003) for information on 
the differences between docking and berthing. For one, 
future research will concern computing optimal berthing 
maneuvers, but also to derive low cost feedback controllers 
for both docking and berthing situations. 



tion of an interior-point filter line-search algorithm for 
large-scale nonlinear programming. Mathematical Pro- 
gramming, 106, 25-57. doi:10.1007/sl0107-004-0559-y. 
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